###########################################
#       Assignment: WGCNA
#       Date: Nov 17,2019   
#       Author: Shawn Wang
###########################################

##=====01.准备========

setwd("/Volumes/FileManage/0660Drive/02.heterosis_project/08133BtGroup/03.WGCNA/04.newVersion/")
options(stringsAsFactors = F)
# if (!requireNamespace("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("edgeR")
# BiocManager::install("dplyr")
# BiocManager::install("WGCNA")
# BiocManager::install("reshape2")
# BiocManager::install("tidyverse")
# BiocManager::install("FactoMineR")
# BiocManager::install("factoextra")
# BiocManager::install("pheatmap")
library(edgeR)
library(WGCNA)
library(dplyr)
library(reshape2)
library(stringr)
library(tidyverse)
library(FactoMineR)
library(factoextra)
library(pheatmap)
source("/Volumes/FileManage/0660Drive/02.heterosis_project/09.code/11.01.PCAandHeatmap.R")
source("/Volumes/FileManage/0660Drive/02.heterosis_project/09.code/11.02.WGCNA.SFT.R")
source("/Volumes/FileManage/0660Drive/02.heterosis_project/09.code/11.03.WGCNA.module.R")
source("/Volumes/FileManage/0660Drive/02.heterosis_project/09.code/11.04.WGCNA.moduleTrait.R")
source("/Volumes/FileManage/0660Drive/02.heterosis_project/09.code/11.05.WGCNA.HubGene.R")
# 原始count数据
rawCount <- read.table("/Volumes/FileManage/0660Drive/02.heterosis_project/08133BtGroup/03.WGCNA/02.Rawdata/mRNA_readcount.xls",
                       header = T,
                       sep = "\t")

# head(rawCount)
rownames(rawCount) <- rawCount$transcript_id
rawCount <- rawCount[,-1]
bt.raw <- data.frame(row.names = rownames(rawCount),
                     select(rawCount, contains("Bt_")),
                     select(rawCount, contains("J_")),
                     select(rawCount, contains("SY_")))
bt.raw <- select(bt.raw, -contains("Z12"))

# 按照组织分成3份
leaf <- select(bt.raw, contains("_L"))
ovule <- select(bt.raw, contains("_O"))
fiber <- select(bt.raw, contains("_F"))

##======02.方法======================
# 参数设置
type = "unsigned"
corType = "pearson"
corFnc = ifelse(corType=="pearson", cor, bicor)
maxPOutliers = ifelse(corType=="pearson",1,0.05)
robustY = ifelse(corType=="pearson",T,F)
# 组织
tissue = ovule
# 名字
Title = "ovule"

##======02 样品PCA和Heatmap============
WGCNA.PCAandHeatmap(tissue = tissue,
                    Title = Title)

##======02 软阈值SFT计算===============
# 表达
datExpr = y
WGCNA.SFT(Title = Title,
          datExpr = datExpr)

##======02 Module计算===============
# 无向网络在power小于15或有向网络power小于30内，没有一个power值可以使
# 无标度网络图谱结构R^2达到0，平均连接度较高如在100以上，可能是由于
# 部分样品与其他样品差别太大。这可能由批次效应、样品异质性或实验条件对
# 表达影响太大等造成。可以通过绘制样品聚类查看分组信息和有无异常样品。
# 如果这确实是由有意义的生物变化引起的，也可以使用下面的经验power值。
if (is.na(power)){
power = ifelse(nSamples<20, ifelse(type == "unsigned", 9, 18),
               ifelse(nSamples<30, ifelse(type == "unsigned", 8, 16),
                      ifelse(nSamples<40, ifelse(type == "unsigned", 7, 14),
                             ifelse(type == "unsigned", 6, 12))
               )
)
}

WGCNA.oneStepNetWork(Title = Title)

##=======表型鉴定========
## 分开做4133Bt组合
rownames(datExpr)
phenotype = data.frame(row.names = rownames(datExpr),
                       high = rep(c(0,1,0,1,0),times = c(3,3,3,3,3)),
                       median = rep(c(1,0), times = c(3,12)),
                       low = rep(c(0,1,0,1), times = c(6,3,3,3)))
## 总体做产量
# unique(gsub("_.*","",rownames(datExpr)))
# phenotype = data.frame(row.names = rownames(datExpr),
#                        High = rep(c(1,1,1,1,1,0,0,1,1,0,0,1,1,0,1,1,0),each = 3),
#                        Low = rep(c(0,0,0,0,0,0,1,0,0,1,1,0,0,1,0,0,1),each = 3))
WGCNA.ModuleTrait(Title = Title,
                  phenotype = phenotype)
##=======选择感兴趣的模块分析=========
## design设计
datTraits <- data.frame(row.names = rownames(phenotype),
                        subtype = factor(rep(c("Y","Y_Q","Y","Y_Q","Q"),each = 3),
                                         levels = c("Y","Q","Y_Q")))
design = model.matrix(~0 + datTraits$subtype)
colnames(design) = levels(datTraits$subtype)

## 参数设置
moduleName = "MM.yellow"
phenoName = "Y"
## cor 和性状关联的相关系数
cor = 0
## con 模块内部基因的连接度
con = 00
## 导出hubgene
WGCNA.HubGene(Title = Title,
              cor = cor,
              con = con,
              moduleName = moduleName,
              phenoName = phenoName)
##=======导出cytoscape=========
## 导出cytoscape
## 选择degree的阈值
threshold = 0
## 选择模块名称
module = "blue"

WGCNA.Cytoscape(module = module,
                threshold = threshold)
